In [1]:
# make button for code hiding
from hax.ipython import code_hider
code_hider()
Out[1]:
Introduction¶
Using the same techniques as for the XENON100 multi-event noise analysis, the noise spectra of long baseline data with a window length of 800 µs are combined and analyzed in order to pinpoint possible noise sources.
Dataset¶
| commissioning run | name | events | events used | window length | ZLE | HV |
|---|---|---|---|---|---|---|
| 25 | 160226_1218 |
305 | first 40 | 800 µs | no | on |
Setting up the analysis¶
Before calculating the FFTs, some cross-checks are made to test the dataset quality. More specifically, each event is tested for having a single pulse for each channel which is exactly 800 µs long.
In [1]:
# standard library imports
from os.path import join
# third party imports
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import tables as tb
from tqdm import tqdm
# pax imports
from pax import core
from pax.configuration import load_configuration
from pax.plugins.io.BSON import DecodeZBSON
# magic function to enable inline plots (i.e. plots appear in the notebook itself)
%matplotlib inline
# set global plotting options
plt.rcParams["font.size"] = 20
plt.rcParams["figure.dpi"] = 100
plt.rcParams["figure.figsize"] = (8, 6)
# set basic path constants for easier path changes if needed
BASE_PATH = "/home/cichon/xenon1t/analysis/noise_frequency"
DATA_BASE_PATH = "/archive_lngs100TB/common/daqtest/xenon1t"
# define file path constants for later use
DATA_PATH = join(DATA_BASE_PATH, "160226_1218")
# set up pytables filter for compression
FILTERS = tb.Filters(complib='zlib')
# set up PAX facilities (ignore the warning in the output)
CUSTOM_CONFIG = {"pax": {"plugin_group_names": ["input"],
"input": "Zip.ReadZipped",
"input_name": DATA_PATH}}
pax_config = load_configuration(config_names="XENON100", config_dict=CUSTOM_CONFIG)
custom_pax = core.Processor(config_names="XENON100", config_dict=CUSTOM_CONFIG)
bson_decoder = DecodeZBSON(pax_config["DEFAULT"], custom_pax)
In [2]:
# get first event to extract information about sample length
event_generator = custom_pax.get_events()
event_proxy = next(event_generator)
event = bson_decoder.transform_event(event_proxy)
pulse = event.pulses[0]
# constant definitions
LEN_DATA = len(pulse.raw_data)
LEN_PADDED = 2**(int(np.log2(LEN_DATA)) + 1)
FREQS_PADDED = np.fft.rfftfreq(LEN_PADDED, 1e-8)
In [3]:
# dataset quality check
event_generator = custom_pax.get_events()
pulse_num = []
pulse_len = []
for event_proxy in event_generator:
event = bson_decoder.transform_event(event_proxy)
pulse_num.append(len(event.pulses))
for pulse in event.pulses:
pulse_len.append(len(pulse.raw_data))
pulse_num = np.array(pulse_num)
pulse_len = np.array(pulse_len)
print("Number of events:", len(pulse_num))
print("Maximum number of pulses in an event:", pulse_num.max())
print("Number of events with maximum number of pulses:", len(pulse_num[pulse_num==pulse_num.max()]))
print("Maximum pulse length:", pulse_len.max())
print("Number of pulses with maximum length:", len(pulse_len[pulse_len==pulse_len.max()]))
As it can be seen in the output above, the dataset seems to be fine. Because of the huge amount of FFT coefficients calculated, resulting in a correspondingly huge amount of storage space needed, only the first 40 events are being looked at.
In [4]:
# define pytables column descriptions
class FftCoeffs(tb.IsDescription):
event_nr = tb.UInt32Col()
event_time = tb.UInt64Col()
channel_nr = tb.UInt16Col()
coeffs = tb.ComplexCol(16, shape=len(FREQS_PADDED))
amplitudes = tb.FloatCol(8, shape=len(FREQS_PADDED))
phases = tb.ComplexCol(16, shape=len(FREQS_PADDED))
phases_raw = tb.ComplexCol(16, shape=len(FREQS_PADDED))
phases_custom = tb.ComplexCol(16, shape=len(FREQS_PADDED))
# prepare HDF5 file
FFT_DATA_PATH = join(BASE_PATH, "output/fft_data_xenon1t_160226_1218.h5")
try:
data_file = tb.open_file(FFT_DATA_PATH, mode = "r", filters=FILTERS)
coeffs_table = data_file.get_node("/fft_coeffs")
process_new = False
except (OSError, FileNotFoundError, ValueError, tb.NodeError):
data_file = tb.open_file(FFT_DATA_PATH, mode = "w", filters=FILTERS)
coeffs_table = data_file.create_table("/", "fft_coeffs", FftCoeffs)
process_new = True
In [5]:
# calculate FFT (with average subtraction, Blackman window and zero padding) for all events
# also, determine event start times relative to first event
valid_events = 0
set_continue = False
if process_new:
event_generator = custom_pax.get_events()
for event_nr, event_proxy in enumerate(tqdm(event_generator)):
# only use first 40 valid events due to storage file growing too big otherwise
if valid_events >= 40:
break
event = bson_decoder.transform_event(event_proxy)
# skip event if not all PMTs have a full pulse
if len(event.pulses) != pulse_num.max():
continue
for pulse in event.pulses:
if len(pulse.raw_data) != pulse_len.max():
set_continue = True
break
if set_continue:
set_continue = False
continue
if valid_events == 0:
start_time_first_ev = event.start_time
for channel_nr, pulse in enumerate(event.pulses):
# set FFT metadata
coeffs_table.row["event_nr"] = valid_events
coeffs_table.row["event_time"] = event.start_time - start_time_first_ev
coeffs_table.row["channel_nr"] = channel_nr
# calculate FFT
raw_data_avg = np.average(pulse.raw_data) - pulse.raw_data
raw_data_win = np.blackman(LEN_DATA) * raw_data_avg
coeffs_table.row["coeffs"] = np.fft.rfft(raw_data_win, n=LEN_PADDED) / LEN_PADDED
# calculate amplitudes, phases and corrected phases
# also correct FFT coefficients
coeffs_table.row["amplitudes"] = np.abs(coeffs_table.row["coeffs"])
coeffs_table.row["phases_raw"] = np.exp(complex(0, 1)*np.angle(coeffs_table.row["coeffs"]))
coeffs_table.row["phases"] = (coeffs_table.row["phases_raw"] *
np.exp(complex(0, -1)*2*np.pi*np.mod(FREQS_PADDED*1e-9 *
coeffs_table.row["event_time"], 1)))
coeffs_table.row["phases_custom"] = (coeffs_table.row["phases_raw"] *
np.exp(complex(0, -1)*2*np.pi*np.mod(2.0825*1e-4 *
coeffs_table.row["event_time"], 1)))
coeffs_table.row["coeffs"] *= np.exp(complex(0, -1)*2*np.pi*np.mod(FREQS_PADDED*1e-9 *
coeffs_table.row["event_time"], 1))
coeffs_table.row.append()
# free memory
del raw_data_avg
del raw_data_win
# flush after each event
coeffs_table.flush()
valid_events += 1
# create table indices for HUGE reading speed improvement
coeffs_table.cols.event_nr.create_index()
print("Finished creating event number index.")
coeffs_table.cols.channel_nr.create_index()
print("Finished creating channel number index.")
coeffs_table.flush()
# free memory
del event_generator
del event
In [6]:
# get event and channel indices
events = list(set(coeffs_table.read(field="event_nr")))
events.sort()
channels = list(set(coeffs_table.read(field="channel_nr")))
channels.sort()
# specify data paths
AMPL_CH_MEAN_PATH = join(BASE_PATH, "output/ampl_ch_mean_xenon1t_160226_1218.dat")
COEFFS_AVG_PATH = join(BASE_PATH, "output/coeffs_avg_xenon1t_160226_1218.dat")
COEFFS_STD_PATH = join(BASE_PATH, "output/coeffs_std_xenon1t_160226_1218.dat")
AMPL_MEAN_PATH = join(BASE_PATH, "output/ampl_mean_xenon1t_160226_1218.dat")
AMPL_STD_PATH = join(BASE_PATH, "output/ampl_std_xenon1t_160226_1218.dat")
PHASE_MEAN_PATH = join(BASE_PATH, "output/phase_mean_xenon1t_160226_1218.dat")
PHASE_VAR_PATH = join(BASE_PATH, "output/phase_var_xenon1t_160226_1218.dat")
PHASE_RAW_MEAN_PATH = join(BASE_PATH, "output/phase_raw_mean_xenon1t_160226_1218.dat")
PHASE_RAW_VAR_PATH = join(BASE_PATH, "output/phase_raw_var_xenon1t_160226_1218.dat")
PHASE_CUSTOMCORR_MEAN_PATH = join(BASE_PATH, "output/phase_customcorr_mean_xenon1t_160226_1218.dat")
PHASE_CUSTOMCORR_VAR_PATH = join(BASE_PATH, "output/phase_customcorr_var_xenon1t_160226_1218.dat")
# calculate mean amplitude spectrum over all channels for each event
try:
ampl_ch_mean = np.memmap(AMPL_CH_MEAN_PATH, dtype=np.float32, mode="r", shape=(len(events), len(FREQS_PADDED)))
except (FileNotFoundError, ValueError):
ampl_ch_mean = np.memmap(AMPL_CH_MEAN_PATH, dtype=np.float32, mode="w+", shape=(len(events), len(FREQS_PADDED)))
for event_nr in tqdm(events):
np.nanmean(coeffs_table.read_where("event_nr == {}".format(event_nr), field="amplitudes"), 0,
out=ampl_ch_mean[event_nr])
# calculate mean/standard deviation of coefficient/amplitude/phase spectrum over all events for each channel
try:
coeffs_avg = np.memmap(COEFFS_AVG_PATH, dtype=np.complex128, mode="r", shape=(len(channels), len(FREQS_PADDED)))
coeffs_std = np.memmap(COEFFS_STD_PATH, dtype=np.float64, mode="r", shape=(len(channels), len(FREQS_PADDED)))
ampl_mean = np.memmap(AMPL_MEAN_PATH, dtype=np.float64, mode="r", shape=(len(channels), len(FREQS_PADDED)))
ampl_std = np.memmap(AMPL_STD_PATH, dtype=np.float64, mode="r", shape=(len(channels), len(FREQS_PADDED)))
phases_mean = np.memmap(PHASE_MEAN_PATH, dtype=np.complex128, mode="r", shape=(len(channels), len(FREQS_PADDED)))
phases_var = np.memmap(PHASE_VAR_PATH, dtype=np.float64, mode="r",shape=(len(channels), len(FREQS_PADDED)))
phases_raw_mean = np.memmap(PHASE_RAW_MEAN_PATH, dtype=np.complex128, mode="r",shape=(len(channels),
len(FREQS_PADDED)))
phases_raw_var = np.memmap(PHASE_RAW_VAR_PATH, dtype=np.float64, mode="r", shape=(len(channels), len(FREQS_PADDED)))
phases_customcorr_mean = np.memmap(PHASE_CUSTOMCORR_MEAN_PATH, dtype=np.complex128, mode="r",shape=(len(channels),
len(FREQS_PADDED)))
phases_customcorr_var = np.memmap(PHASE_CUSTOMCORR_VAR_PATH, dtype=np.float64, mode="r", shape=(len(channels),
len(FREQS_PADDED)))
except (FileNotFoundError, ValueError):
coeffs_avg = np.memmap(COEFFS_AVG_PATH, dtype=np.complex128, mode="w+", shape=(len(channels), len(FREQS_PADDED)))
coeffs_std = np.memmap(COEFFS_STD_PATH, dtype=np.float64, mode="w+", shape=(len(channels), len(FREQS_PADDED)))
ampl_mean = np.memmap(AMPL_MEAN_PATH, dtype=np.float64, mode="w+", shape=(len(channels), len(FREQS_PADDED)))
ampl_std = np.memmap(AMPL_STD_PATH, dtype=np.float64, mode="w+", shape=(len(channels), len(FREQS_PADDED)))
phases_mean = np.memmap(PHASE_MEAN_PATH, dtype=np.complex128, mode="w+", shape=(len(channels), len(FREQS_PADDED)))
phases_var = np.memmap(PHASE_VAR_PATH, dtype=np.float64, mode="w+",shape=(len(channels), len(FREQS_PADDED)))
phases_raw_mean = np.memmap(PHASE_RAW_MEAN_PATH, dtype=np.complex128, mode="w+",shape=(len(channels),
len(FREQS_PADDED)))
phases_raw_var = np.memmap(PHASE_RAW_VAR_PATH, dtype=np.float64, mode="w+", shape=(len(channels), len(FREQS_PADDED)))
phases_customcorr_mean = np.memmap(PHASE_CUSTOMCORR_MEAN_PATH, dtype=np.complex128, mode="w+",shape=(len(channels),
len(FREQS_PADDED)))
phases_customcorr_var = np.memmap(PHASE_CUSTOMCORR_VAR_PATH, dtype=np.float64, mode="w+", shape=(len(channels),
len(FREQS_PADDED)))
for channel_nr in tqdm(channels):
data = coeffs_table.read_where("channel_nr == {}".format(channel_nr))
np.nanmean(data["coeffs"], 0, out=coeffs_avg[channel_nr])
np.nanstd(data["coeffs"], 0, out=coeffs_std[channel_nr])
np.nanmean(data["amplitudes"], 0, out=ampl_mean[channel_nr])
np.nanstd(data["amplitudes"], 0, out=ampl_std[channel_nr])
np.nanmean(data["phases"], 0, out=phases_mean[channel_nr])
np.nanvar(data["phases"], 0, out=phases_var[channel_nr])
np.nanmean(data["phases_raw"], 0, out=phases_raw_mean[channel_nr])
np.nanvar(data["phases_raw"], 0, out=phases_raw_var[channel_nr])
np.nanmean(data["phases_custom"], 0, out=phases_customcorr_mean[channel_nr])
np.nanvar(data["phases_custom"], 0, out=phases_customcorr_var[channel_nr])
# translate phase factor means / variances into equivalent phases (in periods)
periods_mean = np.mod(np.angle(phases_mean)/(2*np.pi), 1)
periods_std = np.arccos(1-phases_var/2) / (2*np.pi)
periods_raw_mean = np.mod(np.angle(phases_raw_mean)/(2*np.pi), 1)
periods_raw_std = np.arccos(1-phases_raw_var/2) / (2*np.pi)
periods_customcorr_mean = np.mod(np.angle(phases_customcorr_mean)/(2*np.pi), 1)
periods_customcorr_std = np.arccos(1-phases_customcorr_var/2) / (2*np.pi)
First look at amplitude spectrum¶
As seen in Jelle's first look at long-baseline XENON1T data, two noise components which are already known from the XENON100 multi-event noise analysis return (Fig. 01):
- 25 kHz (sinusoidal noise)
- multiples of approx. 207-208 kHz (delta noise, very likely from HV module)
In addition, some other features can be made out:
- broad peak at approx. 700 kHz
- even broader peak between 2.5 MHz and 3 MHz
- faint peak at about 600 kHz (will be important later on)
Figure 01: Amplitude spectrum (averaged over all channels and events, different zoom levels)¶
In [7]:
# calculate mean over all events
ampl_ch_ev_mean = np.mean(ampl_ch_mean, 0)
# make spectrum plots
fig_avg, axes_avg = plt.subplots(figsize=(16, 18), nrows=3, ncols=1)
avg_plot_1 = axes_avg[0].plot(FREQS_PADDED/1e6, ampl_ch_ev_mean)
axes_avg[0].set_xlim(0, 50)
axes_avg[0].set_ylim(1e-3, 1)
axes_avg[0].set_yscale("log")
axes_avg[0].set_title("Event 0 amplitude spectrum (amplitudes averaged over all channels)", y=1.05)
axes_avg[0].set_xlabel("frequency [MHz]")
axes_avg[0].set_ylabel("Amplitude [ADC units]")
axes_avg[0].minorticks_on()
avg_plot_2 = axes_avg[1].plot(FREQS_PADDED/1e6, ampl_ch_ev_mean)
axes_avg[1].set_xlim(0, 10)
axes_avg[0].set_ylim(1e-3, 1)
axes_avg[1].set_yscale("log")
axes_avg[1].set_title("Event 0 amplitude spectrum (amplitudes averaged over all channels)", y=1.05)
axes_avg[1].set_xlabel("frequency [MHz]")
axes_avg[1].set_ylabel("Amplitude [ADC units]")
axes_avg[1].minorticks_on()
avg_plot_3 = axes_avg[2].plot(FREQS_PADDED/1e6, ampl_ch_ev_mean)
axes_avg[2].set_xlim(0, 4)
axes_avg[2].set_ylim(1e-3, 1)
axes_avg[2].set_yscale("log")
axes_avg[2].set_title("Event 0 amplitude spectrum (amplitudes averaged over all channels)", y=1.05)
axes_avg[2].set_xlabel("frequency [MHz]")
axes_avg[2].set_ylabel("Amplitude [ADC units]")
axes_avg[2].minorticks_on()
fig_avg.tight_layout()
Figure 02: Evolution of amplitude spectrum (averaged over all channels, different zoom levels)¶
In [8]:
# set up grid points for plotting
SPECTR_X1, SPECTR_Y1 = np.meshgrid(FREQS_PADDED/1e6, np.arange(0, len(channels)))
SPECTR_X2, SPECTR_Y2 = np.meshgrid(FREQS_PADDED/1e6, np.arange(0, len(events)))
# make evolution plots
fig_avg_evol, axes_avg_evol = plt.subplots(figsize=(16, 18), nrows=3, ncols=1)
avg_evol_plot_1 = axes_avg_evol[0].pcolormesh(SPECTR_X2, SPECTR_Y2, ampl_ch_mean, cmap=mpl.cm.plasma,
norm=mpl.colors.LogNorm(), vmin=1e-3, vmax=3e-1)
axes_avg_evol[0].set_xlim(0, 50)
axes_avg_evol[0].set_ylim(0, 39)
avg_evol_1_cbar = fig_avg_evol.colorbar(avg_evol_plot_1, ax=axes_avg_evol[0])
axes_avg_evol[0].set_title("Evolution of amplitude spectrum (amplitudes averaged over all channels)", y=1.05)
axes_avg_evol[0].set_xlabel("frequency [MHz]")
axes_avg_evol[0].set_ylabel("event number")
axes_avg_evol[0].minorticks_on()
avg_evol_1_cbar.set_label("amplitude [ADC units]", labelpad=22, rotation=270)
avg_evol_1_cbar.ax.minorticks_on()
avg_evol_plot_2 = axes_avg_evol[1].pcolormesh(SPECTR_X2, SPECTR_Y2, ampl_ch_mean, cmap=mpl.cm.plasma,
norm=mpl.colors.LogNorm(), vmin=1e-3, vmax=3e-1)
axes_avg_evol[1].set_xlim(0, 4)
axes_avg_evol[1].set_ylim(0, 39)
avg_evol_2_cbar = fig_avg_evol.colorbar(avg_evol_plot_2, ax=axes_avg_evol[1])
axes_avg_evol[1].set_title("Evolution of amplitude spectrum (amplitudes averaged over all channels)", y=1.05)
axes_avg_evol[1].set_xlabel("frequency [MHz]")
axes_avg_evol[1].set_ylabel("event number")
axes_avg_evol[1].minorticks_on()
avg_evol_2_cbar.set_label("amplitude [ADC units]", labelpad=22, rotation=270)
avg_evol_2_cbar.ax.minorticks_on()
avg_evol_plot_3 = axes_avg_evol[2].pcolormesh(SPECTR_X2, SPECTR_Y2, ampl_ch_mean, cmap=mpl.cm.plasma,
norm=mpl.colors.LogNorm(), vmin=1e-3, vmax=3e-1)
axes_avg_evol[2].set_xlim(0, 1)
axes_avg_evol[2].set_ylim(0, 39)
avg_evol_3_cbar = fig_avg_evol.colorbar(avg_evol_plot_3, ax=axes_avg_evol[2])
axes_avg_evol[2].set_title("Evolution of amplitude spectrum (amplitudes averaged over all channels)", y=1.05)
axes_avg_evol[2].set_xlabel("frequency [MHz]")
axes_avg_evol[2].set_ylabel("event number")
axes_avg_evol[2].minorticks_on()
avg_evol_3_cbar.set_label("amplitude [ADC units]", labelpad=22, rotation=270)
avg_evol_3_cbar.ax.minorticks_on()
fig_avg_evol.tight_layout()
Noise spectrum averages / standard deviations (over all events)¶
In contrast to XENON100 data, there seem to be no noise components which are limited to certain channel ranges only (Fig. 03 and 04, upper halves). In addition, no phase correlation can be made out when combining events (Fig. 03 and 04, lower halves).
Figure 03: Average amplitude and phase spectra¶
In [9]:
# make mean value plots
fig_coeffs, axes_coeffs = plt.subplots(figsize=(16, 32), nrows=5, ncols=1)
coeffs_amplplot_meanlast = axes_coeffs[0].pcolormesh(SPECTR_X1, SPECTR_Y1, ampl_mean, cmap=mpl.cm.plasma,
norm=mpl.colors.LogNorm(), vmin=1e-3, vmax=1)
coeffs_amplplot_meanlast_cbar = fig_coeffs.colorbar(coeffs_amplplot_meanlast, ax=axes_coeffs[0])
axes_coeffs[0].set_xlim(0, 4)
axes_coeffs[0].set_ylim(0, 253)
axes_coeffs[0].set_title("Amplitude spectrum (mean over all amplitudes)")
axes_coeffs[0].set_xlabel("frequency [MHz]")
axes_coeffs[0].set_ylabel("channel number")
axes_coeffs[0].minorticks_on()
coeffs_amplplot_meanlast_cbar.set_label("amplitude [ADC units]", labelpad=20, rotation=270)
coeffs_amplplot_meanlast_cbar.ax.minorticks_on()
coeffs_amplplot_meanfirst = axes_coeffs[1].pcolormesh(SPECTR_X1, SPECTR_Y1, np.abs(coeffs_avg),
cmap=mpl.cm.plasma, norm=mpl.colors.LogNorm(), vmin=1e-7, vmax=1e-1)
coeffs_amplplot_meanfirst_cbar = fig_coeffs.colorbar(coeffs_amplplot_meanfirst, ax=axes_coeffs[1])
axes_coeffs[1].set_xlim(0, 4)
axes_coeffs[1].set_ylim(0, 253)
axes_coeffs[1].set_title("Amplitude spectrum (mean over all FFT coefficients)")
axes_coeffs[1].set_xlabel("frequency [MHz]")
axes_coeffs[1].set_ylabel("channel number")
axes_coeffs[1].minorticks_on()
coeffs_amplplot_meanfirst_cbar.set_label("amplitude [ADC units]", labelpad=20, rotation=270)
coeffs_amplplot_meanfirst_cbar.ax.minorticks_on()
coeffs_phaseplot_meanlast = axes_coeffs[2].pcolormesh(SPECTR_X1, SPECTR_Y1, periods_mean, cmap=mpl.cm.hsv, vmin=0,
vmax=1)
coeffs_phaseplot_meanlast_cbar = fig_coeffs.colorbar(coeffs_phaseplot_meanlast, ax=axes_coeffs[2])
axes_coeffs[2].set_xlim(0, 4)
axes_coeffs[2].set_ylim(0, 253)
axes_coeffs[2].set_title("Phase spectrum (mean over all phase factors)")
axes_coeffs[2].set_xlabel("frequency [MHz]")
axes_coeffs[2].set_ylabel("channel number")
axes_coeffs[2].minorticks_on()
coeffs_phaseplot_meanlast_cbar.set_label("phase [periods]", labelpad=20, rotation=270)
coeffs_phaseplot_meanlast_cbar.ax.minorticks_on()
coeffs_phaseplot_custom = axes_coeffs[3].pcolormesh(SPECTR_X1, SPECTR_Y1, periods_customcorr_mean, cmap=mpl.cm.hsv, vmin=0,
vmax=1)
coeffs_phaseplot_custom_cbar = fig_coeffs.colorbar(coeffs_phaseplot_custom, ax=axes_coeffs[3])
axes_coeffs[3].set_xlim(0, 4)
axes_coeffs[3].set_ylim(0, 253)
axes_coeffs[3].set_title("Phase spectrum (mean over all phase factors; 207 kHz corr.)")
axes_coeffs[3].set_xlabel("frequency [MHz]")
axes_coeffs[3].set_ylabel("channel number")
axes_coeffs[3].minorticks_on()
coeffs_phaseplot_custom_cbar.set_label("phase [periods]", labelpad=20, rotation=270)
coeffs_phaseplot_custom_cbar.ax.minorticks_on()
coeffs_phaseplot_meanfirst = axes_coeffs[4].pcolormesh(SPECTR_X1, SPECTR_Y1,
np.mod(np.angle(coeffs_avg)/(2*np.pi), 1), cmap=mpl.cm.hsv,
vmin=0, vmax=1)
coeffs_phaseplot_meanfirst_cbar = fig_coeffs.colorbar(coeffs_phaseplot_meanfirst, ax=axes_coeffs[4])
axes_coeffs[4].set_xlim(0, 4)
axes_coeffs[4].set_ylim(0, 253)
axes_coeffs[4].set_title("Phase spectrum (mean over all FFT coefficients)")
axes_coeffs[4].set_xlabel("frequency [MHz]")
axes_coeffs[4].set_ylabel("channel number")
axes_coeffs[4].minorticks_on()
coeffs_phaseplot_meanfirst_cbar.set_label("phase [period equiv.]", labelpad=20, rotation=270)
coeffs_phaseplot_meanfirst_cbar.ax.minorticks_on()
fig_coeffs.tight_layout()
Figure 04: Average amplitude and phase spectra (zoomed in)¶
In [10]:
axes_coeffs[0].set_xlim(0, 1)
axes_coeffs[1].set_xlim(0, 1)
axes_coeffs[2].set_xlim(0, 1)
axes_coeffs[3].set_xlim(0, 1)
axes_coeffs[4].set_xlim(0, 1)
fig_coeffs
Out[10]:
The standard deviation plots (Fig. 05 and 06) show, that the 25 kHz noise component cannot assumed to be constant. While the multiples of the delta noise line have a lower standard deviation compared to the rest, there seems to be a certain influence of dark pulses (horizontal lines). Except for the uncorrected phase of the 600 kHz line (Fig. 06, bottom), which can thus be assumed to be readout-related, there are no hints of any phase varying less compared to areas without a line.
Figure 05: Standard deviations of FFT coefficient parameters¶
In [11]:
# make standard deviation plots
fig_coeffs_std, axes_coeffs_std = plt.subplots(figsize=(16, 32), nrows=5, ncols=1)
coeffs_coeffplot_std = axes_coeffs_std[0].pcolormesh(SPECTR_X1, SPECTR_Y1, coeffs_std/ampl_mean, cmap=mpl.cm.plasma,
norm=mpl.colors.LogNorm(), vmin=1, vmax=10)
coeffs_coeffplot_std_cbar = fig_coeffs_std.colorbar(coeffs_coeffplot_std, ax=axes_coeffs_std[0])
axes_coeffs_std[0].set_xlim(0, 4)
axes_coeffs_std[0].set_ylim(0, 253)
axes_coeffs_std[0].set_title("FFT coefficients (standard deviation)")
axes_coeffs_std[0].set_xlabel("frequency [MHz]")
axes_coeffs_std[0].set_ylabel("channel number")
axes_coeffs_std[0].minorticks_on()
coeffs_coeffplot_std_cbar.set_label("$\sigma / \mu_{\mathsf{ampl.}}$", labelpad=23, rotation=270)
coeffs_coeffplot_std_cbar.ax.minorticks_on()
coeffs_amplplot_std = axes_coeffs_std[1].pcolormesh(SPECTR_X1, SPECTR_Y1, ampl_std/ampl_mean, cmap=mpl.cm.plasma,
norm=mpl.colors.LogNorm(), vmin=1e-1, vmax=10)
coeffs_amplplot_std_cbar = fig_coeffs_std.colorbar(coeffs_amplplot_std, ax=axes_coeffs_std[1])
axes_coeffs_std[1].set_xlim(0, 4)
axes_coeffs_std[1].set_ylim(0, 253)
axes_coeffs_std[1].set_title("Amplitudes (standard deviation)")
axes_coeffs_std[1].set_xlabel("frequency [MHz]")
axes_coeffs_std[1].set_ylabel("channel number")
axes_coeffs_std[1].minorticks_on()
coeffs_amplplot_std_cbar.set_label("$\sigma / \mu$", labelpad=23, rotation=270)
coeffs_amplplot_std_cbar.ax.minorticks_on()
coeffs_phaseplot_std = axes_coeffs_std[2].pcolormesh(SPECTR_X1, SPECTR_Y1, periods_std, cmap=mpl.cm.plasma, vmin=0.155)
coeffs_phaseplot_std_cbar = fig_coeffs_std.colorbar(coeffs_phaseplot_std, ax=axes_coeffs_std[2])
axes_coeffs_std[2].set_xlim(0, 4)
axes_coeffs_std[2].set_ylim(0, 253)
axes_coeffs_std[2].set_title("Phases (standard deviation)")
axes_coeffs_std[2].set_xlabel("frequency [MHz]")
axes_coeffs_std[2].set_ylabel("channel number")
axes_coeffs_std[2].minorticks_on()
coeffs_phaseplot_std_cbar.set_label("$\sigma$ [period equiv.]", labelpad=23, rotation=270)
coeffs_phaseplot_std_cbar.ax.minorticks_on()
coeffs_phaseplot_custom_std = axes_coeffs_std[3].pcolormesh(SPECTR_X1, SPECTR_Y1, periods_customcorr_std, cmap=mpl.cm.plasma,
vmin=0.155)
coeffs_phaseplot_custom_std_cbar = fig_coeffs_std.colorbar(coeffs_phaseplot_custom_std, ax=axes_coeffs_std[3])
axes_coeffs_std[3].set_xlim(0, 4)
axes_coeffs_std[3].set_ylim(0, 253)
axes_coeffs_std[3].set_title("Phases (standard deviation; 207 kHz corr.)")
axes_coeffs_std[3].set_xlabel("frequency [MHz]")
axes_coeffs_std[3].set_ylabel("channel number")
axes_coeffs_std[3].minorticks_on()
coeffs_phaseplot_custom_std_cbar.set_label("$\sigma$ [period equiv.]", labelpad=23, rotation=270)
coeffs_phaseplot_custom_std_cbar.ax.minorticks_on()
coeffs_phaseplot_raw_std = axes_coeffs_std[4].pcolormesh(SPECTR_X1, SPECTR_Y1, periods_raw_std, cmap=mpl.cm.plasma,
vmin=0.155)
coeffs_phaseplot_raw_std_cbar = fig_coeffs_std.colorbar(coeffs_phaseplot_raw_std, ax=axes_coeffs_std[4])
axes_coeffs_std[4].set_xlim(0, 4)
axes_coeffs_std[4].set_ylim(0, 253)
axes_coeffs_std[4].set_title("Uncorrected phases (standard deviation)")
axes_coeffs_std[4].set_xlabel("frequency [MHz]")
axes_coeffs_std[4].set_ylabel("channel number")
axes_coeffs_std[4].minorticks_on()
coeffs_phaseplot_raw_std_cbar.set_label("$\sigma$ [period equiv.]", labelpad=23, rotation=270)
coeffs_phaseplot_raw_std_cbar.ax.minorticks_on()
fig_coeffs_std.tight_layout()
Figure 06: Standard deviations of FFT coefficient parameters (zoomed in)¶
In [12]:
axes_coeffs_std[0].set_xlim(0, 1)
axes_coeffs_std[1].set_xlim(0, 1)
axes_coeffs_std[2].set_xlim(0, 1)
axes_coeffs_std[3].set_xlim(0, 1)
axes_coeffs_std[4].set_xlim(0, 1)
fig_coeffs_std
Out[12]:
In [ ]: